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Abstract 

To address the impact of electron correlations in the linear and non-linear response regimes of interacting 
many-electron systems exposed to time-dependent external fields, we study one-dimensional (ID) systems 
where the interacting problem is solved exactly by exploiting the mapping of the ID iV-electron problem 
onto an A^-dimensional single electron problem. We analyze the performance of the recently derived ID local 
density approximation as well as the exact-exchange orbital functional for those systems. We show that the 
interaction with an external resonant laser field shows Rabi oscillations which are detuned due to the lack 
of memory in adiabatic approximations. To investigate situations where static correlations play a role, we 
consider the time-evolution of the natural occupation numbers associated to the reduced one-body density 
matrix. Those studies shed light on the non-locality and time-dependence of the exchange and correlation 
functionals in time-dependent density and density-matrix functional theories. 



1. Introduction 

Since its invention in 1984 timc-dcpcndcnt density- 
functional theory (TDDFT) has become one of the 
major tools for describing time-dependent phenom- 
ena of electronic systems [1, 2]. Despite its success, 
several important questions remain open. A promi- 
nent example are double excitations [3] , which cannot 
be described with adiabatic approximations to the 
exchange-correlation (xc) kernel [4]. Other examples 
include the description of memory [5] , charge-transfer 
excitations [6], Rabi oscillations [7], and population 
control [8, 9]. Also, the construction of function- 
als for certain observables can be problematic, like 
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e.g. double-ionization in strong laser fields where one 
strategy rests on expressing the pair-correlation func- 
tion as a functional of the time-dependent density 
[10]. 

In many cases, there is little knowledge about how 
the dynamics of the many-body system interacting 
with an arbitrary external time-dependent field is 
mapped onto the non-interacting (time-dependent) 
Kohn-Sham system. Here, one-dimensional systems 
can provide insight since these systems can be exactly 
diagonalized and subsequently propagated in time for 
a small number of electrons. We provide insight into 
the limitations of adiabatic functionals, especially for 
describing non-linear electron dynamics exemplified 
by the case of Rabi oscillations. 

This article is organized as follows, we first high- 
light the exact mapping of a many-electron system 
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onto an TV-dimensional one-electron problem and the 
selection of proper fermionic solutions. Then, we 
discuss the recently developed onc-dimensional local 
density approximation (LDA) and its performance 
for calculating linear and non-linear response. We 
use the LDA as well as exact exchange (EXX) to in- 
vestigate the description of double excitations and 
Rabi oscillations with adiabatic approximations. We 
then change from TDDFT to reduced density-matrix 
functional theory, where we discuss under which con- 
ditions adiabatic approximations provide a valid de- 
scription. We conclude the paper with a short sum- 
mary and perspectives. 



2. One-dimensional model systems 



The Hamiltonian for N electrons moving in a gen- 
eral, possibly time dependent, external potential Woxt 
in one spatial dimension reads 



H = 



N 

E 



2dx'j 



j,k=l 



(1) 

where Wmt describes the electron-electron interaction 
(atomic units e = m = U = 1 are used through- 
out this paper). In one spatial dimension the singu- 
larity of the ordinary Coulomb interaction prevents 
electrons from passing the position of the singularity, 
both in the attractive and repulsive case. In order 
to avoid this unphysical behavior of the full Coulomb 
interaction we employ the so called soft-Coulomb in- 
teraction 



Usoft-c(2^i,a;2) 
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+ (a-i - X2)2 



(2) 



instead [11]. Here, qi and q2 describe the charges of 
the particles while a is the usual softening parameter. 
We use a = 1 for all our calculations. Mathemati- 
cally, it is straightforward to show that the Hamilto- 
nian (1) is equivalent to a Hamiltonian for a single 
particle in N dimensions moving in an external po- 



tential 



WNdim(a;i...X7v) 



N 



^Vcxt{Xj) + 

3 = 1 



- y 

j,k=l 
j^k 



Vint{Xj,Xk) 



(3) 

consisting of all the contributions from Wext and Wint. 
The corresponding Schrodinger equation can, hence, 
be solved by any code which is able to treat non- 
interacting particles in the correct number of dimen- 
sions in an arbitrary external potential. 

Due to the Hamiltonian being symmetric under 
particle interchange, Xj O Xk, the solutions of the 
Schrodinger equation can be classified according to 
irreducible representations of the permutation group. 
For the simplest case of two interacting electrons both 
the symmetric and antisymmetric solutions arc valid 
corresponding to the singlet and triplet spin config- 
urations, respectively. For more than two electrons 
one needs to separately ensure that the spatial wave 
function is a solution to the TV-electron problem. For 
example, a totally symmetric spatial wave function 
is a correct solution for a single particle in N dimen- 
sions, however, for N > 2 there is no corresponding 
spin function such that the total wave function has 
the required antisymmetry to be a solution of the 
N particle problem in ID. We solve this problem by 
symmetrizing the solutions according to all possible 
fermionic Young diagrams for the given particle num- 
ber [12]. Fig. 1 shows all possible standard Young 
diagrams for the spatial part of the wave function 
for two and three electrons. As the spin of the elec- 
tron is 1/2, the Young diagrams for the spin part 
can maximally have two rows, one for each spin di- 
rection. The Young diagrams for the spatial part of 
the wave function are then given as the transpose of 
the respective spin diagram and, hence, have at most 
two columns. For two electrons there exist two di- 
agrams corresponding to the singlet (Fig. 1 a) and 
triplet (Fig. 1 b) configurations. For three electrons, 
there exist two diagrams with two electrons in one 
spin channel and the remaining electron in the other 
channel (Fig. 1 c and d) and one diagram with all 
electrons having the same spin (Fig. 1 e). 

In practice, we solve the Schrodinger equation in 
N dimensions and then symmetrize each solution ac- 
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Figure 1: Possible standard Young diagrams for the spatial part of the wave function for two [figures a) and b)] and three 
[figures c) to e)] electrons. There are maximally two columns in each diagram, one for each spin direction. Figures a) and b) 
correspond to the two electron singlet and triplet, respectively. For diagram c) the wave function is symmetrized for particles 
1 and 2 and antisymmetrized for 1 and 3, while for diagram d) the symmetrization is with respect to particles 1 and 3 and 
the antisymmetrization with respect to 1 and 2. For diagram e) the wave function is antisymmetrized with respect to the 
interchange of any two particles. 



cording to the Young diagrams for the given parti- 
cle number. If none of the Young diagrams yields a 
non- vanishing solution after symmetrization the state 
does not describe a solution for spin-1/2 particles 
and is discarded. If a state yields a non-vanishing 
contribution for a given diagram the appropriately 
symmetrized state is normalized and used in further 
calculations. 

The solution of higher dimensional problems within 
these symmetry restrictions has been implemented in 
the octopus computer program [13, 14]. The lowest 
energy solution is found to be purely symmetric and 
is, therefore, for N > 2, discarded. With increas- 
ing number of electrons we also observe an increasing 
number of states which do not satisfy the fermionic 
symmetry requirements. 



3. Local density approximation 



The local density approximation for electrons in- 
teracting in one spatial dimension is derived from 
quantum Monte-Carlo calculations for a ID homo- 
geneous electron gas where the electrons interact via 
the soft-Coulomb interaction in Eq. (2) [15]. The cor- 
relation energy is parametrized in terms of rg and the 
spin polarization ( = {N^ — Ni)/N in the form 

e,(r„C) = ec(r.,C = 0) (4) 
-^C' [ec(r„C = l)-ec(r,s,C = 0)] 



with 

e,xrs,(, - 2A + Brs+Cr^,+Drl 

X ln(l + ars + /30 (5) 

which proves to be very accurate in the parameteri- 
zation for ID systems for different long-range inter- 
actions [16, 17, 15]. Note, that the above energy is 
given in Hartree units. To obtain a priori the exact 
high-density result known from the random-phase ap- 
proximation, i.e. 

ee(r, ->0,C = 0) = -^rl (6) 

e,(r, ->0,C = 1) = -^^^1 (7) 

to leading order in r^, we fix the ratio a/ A to he equal 
to 8/(7r4a2) and l/in^a^) for ( = and C = 1, re- 
spectively. In both cases m is limited to values larger 
than 1. As a result, the number of independent pa- 
rameters in Eq. (5) is reduced to 7. In addition, for 
a = 1 the denominator can be simplified by setting 
B = 0. However, for smaller values of the softening 
parameter the linear term in the denominator is im- 
portant for achieving agreement with the quantum 
Monte-Carlo results. The optimal values of the pa- 
rameters are given in Tab. 1. For more details on 
the ID QMC methodology and the parameterization 
procedure we refer to Refs. [16, 17]. 

We have implemented the ID LDA for a = 1 
in both unpolarized and polarized versions in the 
octopus program [13, 14]. 
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Figure 2: Linear (top) and non-linear (bottom) spectra of 
Be^"*" , calculated from the Fourier transformation of the dipole 
moment using a polynomial damping function, comparing the 
exact and the ID LDA calculation. The inset in the bottom 
figure shows a zoom into the region from 2.7 to 3.0 Ha. 





C-0 


C = i 


A 


18.40(29) 


5.24(79) 


B 


0.0 


0.0 


C 


7.501(39) 


1.568(230) 


D 


0.10185(5) 


0.1286(150) 


E 


0.012827(10) 


0.00320(74) 


a 


1.511(24) 


0.0538(82) 


P 


0.258(6) 


1.56(1.31) • 10"^ 


m 


4.424(25) 


2.958(99) 


A 


6.7-10-^ 


3.3 • IQ-^ 



Table 1: Parameterization of the correlation energy of the ID 
homogeneous electron gas for a softening parameter of a = 1, 
spin unpolarized {C, = 0) and fully polarized {C, = 1) cases are 
given. The error in the last digits is given in parenthesis, while 
the average error. A, (in Hartree) in the full density range is 
given in the last row. 

Fig. 2 shows the hnear and non-hnear absorption 
spectra of a ID Be^"'' system, i.e. with an external 
potential of 

containing two electrons. We use the LDA as an 
adiabatic approximation to the exact time-dependent 
exchange-correlation potential. The spectrum is cal- 
culated in linear response to a spatially constant per- 
turbation at t = 0, i.e. we apply an additional exter- 
nal electric field £ in dipole approximation 

v^^{x,t)^xE^5{t) (9) 

which gives an initial momentum to the electrons. 
The time propagations were performed in a box rang- 
ing from -150 to 150 bohr with absorbing boundary 
conditions [13] and a grid spacing of 0.2 bohr for a to- 
tal propagation time of 10'^ a.u. In the linear regime a 
kick of £o = 10""' Ha/bohr was employed which was 
then increased to 0.01 Ha/bohr to obtain the non- 
linear response. The values of the excitation energies 
can be found in Tab. 2. In linear response, we see 
five peaks in the LDA spectrum which compare well 
with the first five excitations in the exact case. As 
expected, the agreement is better for lower lying ex- 
citations and gets worse the closer we get to the con- 
tinuum. As a guide for the eye wc included the KS 
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HOMO energy of the LDA calculation, enoMO = 2.06 
Ha and the exact ionization potential of 2.41 Ha. The 
onset of the continuum itself appears at too low ener- 
gies in the LDA calculation missing two more clearly 
visible peaks in the exact spectrum. In other words, 
the LDA fails to reproduce the proper Rydberg se- 
ries, a behavior well known from 3D calculations. For 
comparison we also included the results from an EXX 
calculation (which for two electrons is adiabatic and 
equal to Hartree-Fock) which shows a slightly better 
agreement than LDA for the first three excitations 
but, more importantly, reproduces the Rydberg series 
due to the correct asymptotic behavior of the corre- 
sponding exchange potential. The quality of the EXX 
results also implies that correlation is of secondary 
importance in the system for a = 1. The non- linear 
spectrum shows the same excitations as the linear 
spectrum and three additional peaks for the exact 
and the EXX calculation and two additional peaks 
in the LDA spectrum. Their energies are also listed 
in Tab. 2. Due to the spatial symmetry of the sys- 
tem all even order responses are zero and the first 
non- vanishing higher-order response is of third order. 
The frequency Oi = 0.28 Ha corresponds to an ex- 
citation from the second to the third excited state, 
where the transition from the ground to the second 
excited state is dipole forbidden and, hence, can only 
be reached in a two-photon process. The other two 
frequencies, = 0.42 Ha and Jla = 0.54 Ha, corre- 
spond to the transitions from first to second and sec- 
ond to fifth excited state, respectively. Again, both 
the EXX and the LDA calculations yield a good de- 
scription of the low lying excitations, only the third 
peak cannot be resolved in the LDA spectrum. 

One feature of the exact spectrum that is missing 
from both the LDA and the EXX spectra is the small 
dip at 2.8 Ha, see inset in Fig. 2. It results from a 
Fano resonance [18, 19], i.e. the decay of an excited 
state into continuum states. It is missing from both 
approximate spectra due to the double-excitation 
character of the involved excited state. Double ex- 
citations in the linear regime can only be described 
in TDDFT if a frequency-dependent xc kernel is em- 
ployed [4]. Any adiabatic approximation, however, 
leads to a frequency independent kernel. Hence, dou- 
ble excitations, as well as any resulting features, are 
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Figure 3: Discrete part of the linear (top) and non-linear (bot- 
tom) spectra of Be+, calculated from the Fourier transforma- 
tion of the dipole moment using an exponential damping func- 
tion, comparing the exact and the ID LDA calculation. The 
exact ionization potential is at 0.83 Ha and the LDA HOMO 
at 0.63 Ha. 



missing from both the ALDA and the AEXX linear 
response spectra. Apart from the well-known short- 
comings of not including double-excitations and not 
giving the correct Rydberg series, the ID ALDA re- 
produces both the linear and the non-linear exact 
spectra quite well. 

Fig. 3 shows the discrete part of the linear and non- 
linear spectra for the Be"*" system, i.e. the external 
potential is given by Eq. (8) and the system contains 
three electrons. The ionization potential for the exact 
calculation is 0.83 Ha which is again underestimated 
by the LDA HOMO energy of 0.62 Ha. For this sys- 
tem the projection onto the Young diagrams becomes 
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CJi UJ2 UJy, Ldi W5 CJg ^^7 


r^i r^2 ^3 


LDA 


1.10 1.74 1.90 1.96 2.00 


0.22 0.40 


EXX 


1.13 1.82 2.08 2.20 2.27 2.30 2.32 


0.26 0.43 0.52 


exact 


1.12 1.81 2.08 2.19 2.26 2.29 2.32 


0.28 0.42 0.54 



Table 2: Excitation energies from linear and non-linear response of the ID Be + atom corresponding to the spectra in Fig. 2. 
Excitations from linear response are denoted as uj while those from the non-linear spectrum are denoted with Q. All numbers 
are given in Hartree. 



important with the lowest energy spatial solution be- 
ing symmetric under exchange of any two variables. 
Therefore, it is not a valid solution for three fermions 
and, hence, discarded. The second lowest energy is 
doubly degenerate with the eigenstates correspond- 
ing to diagrams Ic and Id. One of these states is 
then propagated with a kick strength of Eq = 10~* 
for the linear spectra and Eq = 0.1 for the non- linear 
spectra. The exact linear spectrum shows two transi- 
tions at 0.36 Ha and 0.62 Ha. Again, wc observe that 
LDA underestimates these excitation energies giving 
0.34 Ha and 0.55 Ha, respectively. The non-linear 
spectrum contains two more peaks in the exact spec- 
trum at 0.09 Ha and 0.16 Ha which are, however, 
difficult to resolve. In the LDA spectrum only the 
peak at 0.16 Ha can be resolved. From the exact cal- 
culation of excited states we know that there should 
be several more transitions at very small frequencies 
which are accessible in non-linear response. Those 
arc, however, very close to each other and, hence, 
more difficult to resolve. Attempts to improve the 
spectra in the small frequency region are currently in 
progress. 



4. Double excitations 

In order to investigate double excitations in the 
linear-response spectrum we employ the following ex- 
ternal potential 



, .cosh 



cosh (kx) 



(10) 



for which the one-particle problem can be solved an- 
alytically [12] and the resulting eigenvalues are given 













^ / 
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-8.0 Ha 



Figure 4: External cosh-potential and single particle eigen- 
values. In the non-interacting two-particle ground state both 
electrons occupy the lowest energy level in a singlet configura- 
tion. 



as 




(11) 



for j = 0, 1.... Here, the term in parenthesis needs 
to be positive which restricts the number of bound 
states of the system. In other words, by choosing 
the two parameters vq and k appropriately, one can 
create systems with any number of bound states. For 
our calculations we choose fc = 1 and vq = 10 which 
leads to the four bound single-electron states shown 
in Fig. 4. 

Putting two electrons into our system we calcu- 
late the total energies for non-interacting electrons as 
well as for electrons interacting via the soft-Coulomb 
interaction (2). The results for the first ten states 
are given in Tab. 3. As the energy differences be- 
tween the interacting and non-interacting cases are 
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Symmetric well 


Asymmetric well 




State j 


non- interact. 


interacting 


non- interact. 


interacting 


Spin 





-16.00 


-15.10 


-16.92 


-16.02 


singlet 


1 


-12.50 


-11.75 


-13.21 


-12.45 


singlet 


2 


-12.50 


-11.62 


-13.21 


-12.32 


triplet 


3 


-10.00 


-9.31 


-10.52 


-9.83 


singlet 


4 


-10.00 


-9.30 


-10.52 


-9.81 


triplet 


5 


-9.00 


-8.22 


-9.48 


-8.70 


singlet 


6 


-8.50 


-7.98 


-8.90 


-8.38 


singlet 


7 


-8.50 


-7.97 


-8.90 


-8.38 


triplet 


8 


-8.00 


-7.90 


-8.45 


-8.36 


singlet 


9 


-8.00 


-7.90 


-8.45 


-8.36 


triplet 



Table 3: Energies (in Hartroe) for two-particle states in the symmetric well potential (10) and the asymmetric well Eq. (12) 
without interaction and with soft-Coulomb interaction, Eq. (2). The 5th excited state corresponds to a double excitation. 







Symmetric well 




Asymmetric well 


Excitation 


Symmetry 




non-interact. 


interact. 


non-interact. 


interact. 




odd 




-3.50 


-3.48 


-3.71 


-3.70 




even 




-6.00 


-5.80 


-6.40 


-6.21 


rill 


even 




-7.00 


-6.88 


-7.44 


-7.32 




odd 




-7.50 


-7.13 


-8.02 


-7.64 


^12 


odd 




-9.50 


~ -9.28 


-10.11 


- -9.91 



Table 4: Excitation energies (in Hartree) for two particles in the symmetric well potential (10) and the asymmetric well 
Eq. (12) without interaction and with soft-Coulomb interaction, Eq. (2). On and are double excitations. We also state 
the symmetry of the two-particle excited state for the symmetric well. 
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small we can treat the many-body states as perturbed 
independent-particle states. This treatment is con- 
venient because in the independent-particle picture 
double excitations are well defined: they describe 
transitions in which two electrons get excited with 
the excitation energies given as the sum of two single- 
particle excitations. The excitation energies for the 
different transitions are shown in Tab. 4. We note 
that the two-particle eigenstates of the symmetric 
well (10) can be chosen as eigenstates of the par- 
ity operator and, hence, can be classified as even and 
odd. For odd operators like the dipole operator the 
transitions from the ground state (even) to even two- 
particle excited states have zero oscillator strength. 
Nevertheless, these transitions can be visible beyond 
linear response. In addition, starting from the non- 
interacting ground state, doubly excited states have 
zero weight in the density response function because 
the density operator is a single particle operator. 
Thus, also the odd doubly excited two-particle states 
have zero oscillator strength for non-interacting par- 
ticles. The 5th excited state of the non-interacting 
electrons can clearly be identified as a double exci- 
tation. This transition corresponds to both electrons 
getting promoted to the first excited state e\. As the 
energy differences between the interacting and non- 
interacting cases are small the 5th excited state of 
the interacting system is of double-excitation charac- 
ter as well. Unfortunately, however, this two-particle 
excited state is even under parity and, hence, the 
transition is dipole forbidden. The first double ex- 
citation which is dipole allowed is r2i2 which, in the 
independent-particle picture, corresponds to one elec- 
tron getting promoted to the first excited state e\ 
and the other to the second excited state £2- It has 
an excitation energy of 9.50 Ha in the non-interacting 
system. The first ionization potential of the system, 
however, is 8.00 Ha which implies that the dipole al- 
lowed double excitation lies in the continuum part of 
the spectrum. In addition, r2i2 has zero weight in lin- 
ear response for the non-interacting system because 
the transition matrix element vanishes as the final 
state differs from the initial state in two orbital oc- 
cupations. For the interacting system, however, the 
two-particle spatial singlet wave function is no longer 
given as a product of the lowest energy single-particle 




Energy [a.u.] 

Figure 5: Linear response spectrum for two electrons in a 
cosh potential calculated from the Fourier transformation of 
the dipole moment using a polynomial damping function. For 
interacting electrons we observe an additional transition at 
9.4 Ha that corresponds to f2i2 (see inset). 



orbital. A configuration interaction (CI) expansion 
of this wave function also contains terms which cor- 
respond to single excitations of the non-interacting 
particles. As a result, the double excitation r2i2 be- 
comes accessible in linear response. This can be seen 
in Fig. 5, where we plot the absorption spectrum of 
the two-electron system both for interacting and non- 
interacting electrons. The spectrum was calculated in 
linear response to a spatially constant perturbation 
at t = 0, see Eq. (9). We observe that the interacting 
spectrum shows a small dip at sa 9.4 Ha which is close 
to the energy difference between the ground state and 
the dipole-allowed double excitation described earlier 
(as this excitation lies within the continuum its en- 
ergy cannot be computed directly but an estimate 
can be found from fioi + ilo2)- We can clearly see 
that the transition lies within the continuum, or due 
to the calculation being done in a finite box, within 
the excitations to box states. Due to the nearby ex- 
citations to the continuum the frequency f2i2 of the 
bound transition is shifted slightly [18]. This exci- 
tations appears as a dip rather than a peak in the 
spectrum due to the absorbing boundary conditions 
which were employed in the calculation [19]. 

In order to investigate the double excitations which 



8 



le+04 
le+02 



— non-interacting 

- - soft Coulomb 




4 6 
Energy [a.u.] 



Figure 6: Linear response spectrum for two electrons in 
the modified cosh potential, Eq. (12), calculated from the 
Fourier transformation of the dipole moment using a polyno- 
mial damping function. Notice that for interacting electrons 
the double excitation Sin appears at Rj 7.3 Ha and Q12 shifts 
to 10.0 Ha. 



le+04 
le+02 

■g ie+00 

s 

]i ie-02 

■a le-04 



J= le-06 



le-08 
le-10 



L ' 1 ' 


1 ' 1 




1 ' 1 J 


1 EXX 
= - - exact 
-■ LDA 




















\\ 

1 : 
1 ' 








n 

li 




- , 1 . 


1 , 1 


[i , 


1.1- 







4 6 8 10 

Energy [a.u.] 



Figure 7: Linear response spectrum for two electrons in 
the modified cosh potential, Eq. (12), calculated from the 
Fourier transformation of the dipole moment using a polyno- 
mial damping function. For the DFT (EXX and LDA) spectra 
both double excitations On and f2i2 are missing. 
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Figure 8: Non-linear response spectrum for two electrons using 
EXX in the modified cosh potential, Eq. (12), calculated from 
the Fourier transformation of the dipole moment using a poly- 
nomial damping function. The double excitation appears 
at !5i 7.4 Ha (see arrow). 



are dipole forbidden by symmetry, we break the spa- 
tial symmetry of the system by modifying the exter- 
nal potential to 



wo(l + 0.5x) 
cosh^(fcx) 



(12) 



As the Hamiltonian for this case no longer commutes 
with the parity operator, the eigcnstates do not have 
a specific symmetry any longer. Therefore, the pre- 
viously dipole forbidden transitions now have a finite 
oscillator strength. Also, the modification is small 
enough to leave the ordering of the states intact, i.e. 
the fifth excited state still has double-excitation char- 
acter and the energies are approximately those of the 
symmetric system (see Tab. 3 for details). As we 
can see in Fig. 6 this leads to an additional peak 
slightly above 7.3 Ha which corresponds to the tran- 
sition r^ii from the ground state to the fifth excited 
state. In Fig. 7 we also include DFT linear response 
spectra for the symmetry-broken potential Eq. (12) 
using EXX and the ID LDA functional of section 
3 which are used as adiabatic approximations to the 
time-dependent xc potential. It is, therefore, not sur- 
prising that the double excitation i^n is missing from 
the resulting spectra [4, 20]. However, as we can see 
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in Fig. 8, beyond linear response the double excita- 
tion fill becomes visible even for the EXX functional, 
which is adiabatic, in the symmetry-broken potential. 

5. Rabi oscillations 

Rabi oscillations can occur when a system is ex- 
posed to an external laser field with frequency w, 
which is in resonance to a transition in the sys- 
tem. The system oscillates between two states with 
frequency f2o provided it can approximately be de- 
scribed as a two-level system. This is the case if the 
frequency u is resonant to one specific transition in 
the system and the frequency of the oscillation be- 
tween the states is much smaller than the frequency 
of the applied laser field. As will be discussed later, 
a small detuning of the applied laser from the ex- 
act resonant frequency still leads to Rabi oscillations, 
however, with increased frequency and smaller ampli- 
tude. 

We analyze Rabi oscillations for a ID two-electron 
model (sec section 2) but note that the analysis can 
easily be extended to three dimensions and, using 
the single-pole approximation also to any number of 
electrons [21]. For ease of comparison we choose the 
same model as in [7] with the external potential 

vf^^\x, t) = + x8o sin(^<). (13) 

yx'' + 1 

Eq. (13) describes a ID Helium atom interacting with 
a monochromatic laser field of frequency cj. We de- 
note the eigenstates and eigenvalues of the Hamilto- 
nian (1) with the external potential (13) as tpk and 
Efe, respectively. In order for Rabi's solution to be 
valid, the system under study needs to be an effective 
two-level system reducing the solution space to ipo 
(ground state) and ipi (dipole allowed excited state) 
with eigenenergies eo and ei. The system then has a 
resonance at A = ei — eg. The two-level approxima- 
tion is valid if the two conditions 

— « 1, r^o << w (14) 

are satisfied, where S = lj — A describes the de- 
tuning from the resonance and i^o = c^io^o is the 



Rabi frequency for a resonant laser, with dio = 
{ipol J2j IV'i) the dipole matrix element. In order to 
satisfy the second condition we choose Eq = 0.0125ix) 
in Eq. (13). For the eigenvalues we obtain cq — 
—2.238 Ha and ei — —1.705 Ha which implies a reso- 
nant frequency of A = ei — eo = 0.534 Ha. The static 
dipole matrix element is diQ ~ 1.104. The frequency 
u> of the applied field has been chosen to be close to 
the resonance, i.e. the detuning (5 = a; — A is small. 

For an effective two-level system the time- 
dependent two-electron wave function ip(xi,X2,t) can 
be written as a linear combination of ground and ex- 
cited state, i.e. 

'ip{xi,X2,t) ^ aQ{t)il;o{xi,X2) + ai{t)'ipi{xi,X2) (15) 

with |ao(i)P = no{t) and |ai(t)p = ni{t) being 
the time-dependent level populations of the ground 
and excited states. Normalization of the wave func- 
tions then implies no(t) + ni{t) = 1. As the am- 
plitude So and the frequency uj of the applied field 
£{t) are chosen such that the conditions (14) are ful- 
filled the Hamiltonian (1) can be projected onto a 
2x2 space. The time-dependent Schrodinger equa- 
tion idt\ip{t)) = H\ip{t)) then reduces to a 2 x 2 ma- 
trix equation of the form 

[a,it) )-[ d,o£{t) ei )[ ai(i) J 

(16) 

from which one can derive coupled differential 
equations for the level population ni{t) and the 
dipole moment d{t) = {'i(j{t)\xi + X2\ipit)) = 
2dioRe{aQ{t)ai{t)). For the dipole moment we ob- 
tain 

d{t) ^ 2dio^/noit)ni{t) cos{Lut + d{t)), (17) 

where the Rabi frequency is included in the time 
dependence of the amplitude of the dipole moment 
through the time dependence of no(t) and ni{t), and 
6{t) is the phase difference of ao(i) and ai{t). Fulfill- 
ment of conditions (14) allows for the use of the ro- 
tating wave approximation (RWA) [22] which yields 
the following differential equation for ni (t) 

d^n,it)^-{S' + nl)n,it) + ^nl (18) 
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tions, only studying the level populations gives the 
complete picture. 

For the system specified above a comparison be- 
tween the analytic solution of Eqs. (17), (18) and 
the results of the time-propagation with the octopus 
code [13, 14] shows a perfect agreement, which con- 
firms that the conditions (14) are fulfilled for the cho- 
sen values of Sq and lu. 

The KS Hamiltonian corresponding to (13) is given 

as 

N 

^H"s+T. U^ci^, , t) + ^j^o sin(^t)) , (19) 



Figure 9: Dipole moment (red) and level populations ni{t) 
(solid black lino) and no{t) (dashed black line) from analytic 
solution of (16) using So = 0. 01250; for detuning S = 0.08 Ho 
(0.0006 Ha) (a) and S = 2.2 Qq (0.016 Ha) (b). 



with initial conditions ?t.i(0) = and ni{0) = 0. 
Eq. (18) describes a harmonic oscillator with a restor- 
ing force which increases with increasing detuning d. 
As a result, the frequency of the Rabi oscillations 
increases with increased detuning, while the max- 
imum population of the excited state decreases as 

In Fig. 9 the time-dependent dipole moment d{t) 
and the level populations rio(i) and ni(t) for 6 = 
0.08 rio and S = 2.2 JIq are shown. The effect of the 
detuning manifests itself in an incomplete population 
of the excited state and a consequent decrease in the 
amplitude of the envelope of the dipole moment that 
is proportional to y/riQUi. For small detuning the 
minima and the maxima of ni coincide with minima 
of the envelope, but for larger detuning the dipole 
moment only goes to zero for the minima of ni. In 
Fig. 

reffig:dipolenlLina the detuning is small but non- 
zero, hence, the neck at the minima of hq. The neck 
grows with increasing d and evolves into a maximum 
for Fig. 9b. Thus, the first minimum in Fig. 9b cor- 
responds to one complete cycle and can be identified 
with the second minimum in Fig. 9a. We note that 
looking only at the dipole moment is insufiicient to 
discern between resonant and detuned Rabi oscilla- 



where the static KS Hamiltonian reads 



N 



t{xj) + VYi^c[po]{xj). (20) 



We denote the eigenfunctions of _ff ° with (j)k (x) and 
their eigenvalues as e| . As we arc using a two-electron 
system, a single orbital is doubly occupied in the KS 
system. The time evolution of this orbital follows 
from the KS equation 



idt<t){x,t) = Hs<t>{x,t) 



(21) 



with the initial condition (f){x,t = 0) = 4>o{x). This 
equation is non-linear due to the dependence of the 
Hartree-exchange-correlation potential Uhxc on the 
density, p{x,t) = 2|^(x,t)p. The time-dependent 
dipole moment d(t) is an explicit functional of the 
density, i.e. d{t) = J xp{x, t) dx. The exact KS sys- 
tem reproduces the exact many-body density p{x, t) 
and, hence, the exact dipole moment d{t). However, 
this need not be true for an approximate functional. 
Especially, using adiabatic approximations has been 
shown to have a dramatic effect on the calculated 
density during Rabi oscillations [7]. 

Propagating with the EXX and ALDA (see section 
3) results in the dipole moments shown in Fig. 10. 
The resonant frequencies are calculated from lin- 
ear response which yields cj^^-^^ = 0.476 Ha and 
^Exx _ 549 then apply a laser field in 

analogy to the exact calculation with an amplitude 
of £o = 0.0125ti; using the resonant frequency for each 
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Figure 10: Dipole moment (red) and level populations (solid 
black line) and (dashed black line) for EXX (a) and ALDA 
(b). 



case. Wc observe in the level populations that both 
EXX and ALDA show the characteristic signatures 
of detuned Rabi oscillations despite the fact that the 
applied laser is in resonance with the system. 

To clarify whether Rabi oscillations are well de- 
scribed in the context of adiabatic TDDFT we study 
the Hamiltonian (19) in more detail. The Hartree and 
xc potentials, Whxc[p], render the KS differential equa- 
tion non-linear. More specifically, for an adiabatic 
approximation the potential at time t is a functional 
of the density at this time, i.e. Whxc(i) ~ Whxc[p(i)]- In 
the following, we show that vi^^dt) introduces a de- 
tuning that drives the system out of resonance. We 
again rely on the conditions (14), i.e. describe the 
KS system as an effective two-level system. There- 
fore, the time-dependent orbital (f>{x, t) is given as a 
linear combination of the ground-state KS orbital (po 
and the first excited state orbital (jji 



4>ix,t) = a'o{t)(j)o{x) + al{t)(l3i{x). 



(22) 



Projecting the KS Hamiltonian (19) onto the two- 
level KS space (22) yields the 2x2 matrix 



+ di,£{t)+T.,4t) 

dyn/ 



(23) 



-F^c(t) = (</'o|«^^"WI<^i)- This matrix enters Eq. (16) 
to determine the coefficients ag(t) and a\{t). Com- 
pared to Eq. (16) we notice that each entry contains 
an additional term depending on wj^^c j i-^- both the 
electric field and the KS eigenvalues are modified. 
In order to investigate the consequences of the addi- 
tional terms we use the EXX functional for which rel- 
atively simple analytic expressions for the additional 
matrix elements can be derived. 

For the two-electron case investigated here, the 
Hartree-exchange-correlation potential Vh^'^ {x,t) is 
equal to half the Hartree potential and, hence, given 
as 



EXX 
hxc 



pojx') + 6p{x',t) 
yj{x - x'Y + 



(24) 



where we split the total density p{x,t) into a time- 
independent contribution po{x) = 2|(/)o(x)p and a 
rest Sp{x, t) which can be calculated from Eq. (22) 



5p{x,t) 



2\a{\\\Mx)\^~\M^r) (25) 



The part of Eq. (24) containing po determines fhxc[po] 
while dp results in the additional w^^c • Since 0o s.nd 
(f>i usually have opposite spatial symmetry (if the 
Hamiltonian is symmetric) the first term in Eq. (25) 
is symmetric while the second is antisymmetric. This 
results in the first term only contributing to the di- 
agonal elements of the matrix (23) while the second 
term only contributes to the off-diagonal elements. 
Defining the level populations in the KS system as 
rijit) = |a^(t)p allows us to rewrite the contributions 
to the diagonal terms as 



ef (t) = X.nlit), 



(26) 



with dl„ = (<^i|x|0o), er(t) = (0,l«,rcWl0.>, and 



where, for the EXX approximation, the coefficient Xj 
reads 

//(l^i(^Oi;-|0o(xOni0,(x)P.^,.^, (27) 
' JJ ^{x - x')2 + a2 

For the off-diagonal contribution we recall that 
d''{t) = 2dfoi?e(a5(t)*af (t)), as in the exact case, 
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and rewrite the contribution of v-^^^{t) to the ofF- 
diagonal terms as a coefficient g muhiphed by the 
time-dependent dipole moment 



J^xc{t) = g 



(28) 



^10 



For the two-electron system in the EXX approxima- 
tion g is given as 
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(l)i{x')(t>o{x')(l)Q{x)(t>i{x) 
^J{x - x'Y + a2 



dxdx' . 



(29) 



The coefficient g also enters the calculation of the 
resonant frequencies in linear response. Within the 
single-pole approximation the resonant frequency is 
given as = ^As(As + 2g) which for the EXX 
functional yields ijj^-^-^ = 0.532 Ha. The deviation 
from the frequency calculated from time propagation 
of the Hamiltonian (19) in octopus is of the order 
of 3%, coinciding with the deviation of our system 
from a true two-level system which we estimate from 
(1 — [nf^it) + nf (t))). Using, as in the exact case, the 
RWA we obtain to leading order in A/cjq and .g/wo the 
following equation of motion for the level population 



1 



dtnlit) = - [^\n{{tY + j nl{t) + -17^ (30) 

with = rffo^^o and 7 = A — 2g. Neglecting the 
higher order terms in A/wq and g/ujQ introduces an 
error of about 10% in favor of keeping the equa- 
tion simple while still containing the important phys- 
ical effects. Unlike Eq. (18) which represents a har- 
monic oscillator, Eq. (30) corresponds to an anhar- 
monic quartic oscillator and its solution is given in 
terms of Jacobi elliptic functions [23]. Equivalently, 
Eq. (30) can be integrated numerically. Even though 
the oscillator is no longer harmonic the detuning 
still results in an increase of the restoring force. In 
other words, the adiabatic approximation introduces 
a time-dependent detuning proportional to ^n\{t). 

Using the same ID model system (13) as before 
and diagonalizing the Hamiltonian (19) gives for 
the bare KS eigenvalues e^-^-^ = —0.750 Ha and 
^Exx ^ _o.257 Ha which yields = 0.494 Ha. For 
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Figure 11: Potentials corresponding to the differential equa- 
tions (18) and (30). The dynamical detuning leads to a quartic 
potential which has a similar effect to the potential as a large 
detuning. 



the various matrix elements wc obtain c?f 



0.897, 



g = 0.071, A = -0.125, and 7 = -0.268. As a result 
the detuning is of the order of 5% of the resonant 
frequency, i.e. quite large compared to the detun- 
ing which is necessary to destroy the resonant Rabi 
behavior (see Fig. 9) which explains the results we 
see for the dipole moment and level populations (see 
Fig. 10a). In Fig. 11 we plot the potential correspond- 
ing to the restoring forces in the differential equations 
(18) and (30). As we can see, the dynamical detun- 
ing in the EXX calculation has a similar effect on the 
squeezing of the potential as the large detuning in the 
linear Rabi oscillations. 

The behavior using ALDA is very similar to the one 
for the EXX approximation (see Fig. 10b). The anal- 
ysis, however, is more involved due to the functional 
not being linear in the density. We can conclude that 
any adiabatic functional, even the exact adiabatic one 
[24] , will lead to a detuning in the description of Rabi 
oscillations due to the lack of memory and the fact 
that the exact density changes dramatically during 
the transition. For the exact functional the detuning 
effect is compensated by a dependence on the density 
at previous times. 
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Figure 12: Changes in the dipole moment (green) and the first two natural occupation numbers (red and blue) for the helium 
atom (left) and the hydrogen molecule (right) in linear response to a (5-kick of strength So = 0.05 Ha/bohr. 



6. Reduced density-matrix functional theory 

The use of adiabatic approximations in DFT is 
known to miss certain aspects of tlie interacting sys- 
tem, one of wlrich being tlie double excitations dis- 
cussed in Section 4. Another example are charge- 
transfer excitations which, while appearing in linear 
response, are not described correctly [6]. Recently, it 
has been suggested that using reduced density-matrix 
functional theory (RDMFT) these problems can be 
addressed, even within adiabatic approximations [25]. 

RDMFT uses the one-body density matrix 
(IRDM) 

7i(r,r') =iv//dV2...dV^«'(r,r2...r^)**(r',r2...rAr) 



and are called natural orbitals and occupation num- 
bers, respectively. The exact kinetic energy of a sys- 
tem can be written explicitly in terms of the IRDM 
which presents a major advantage compared to stan- 
dard density functional theory. 

Since in general the interaction energy is known as 
an exact functional of the diagonal 72(r, r'; r, r') of 
the two-body density matrix (2RDM) 



(31) 

as its basic variable. Approximations within this the- 
ory arc usually stated in terms of the eigenfunctions 
(pj{r) and eigenvalues rij of 7i(r, r') which satisfy 



d^r'7i(r,r')<y9j(r') = njipj{r) 



(32) 



72(1-1, 1-2; r'l.r'a) 



N{N - 1) 



(fr3...(frN (33) 



«'(ri, r2, r3...rjv)**(r'i, r'2, r3...rAr), 

RDMFT can be viewed as a way to express the di- 
agonal of the 2RDM as a functional of the IRDM. 
Using this functional, the interaction energy is writ- 
ten typically as a sum of the Hartree energy and the 
exchange-correlation (xc) energy, where only the cor- 
relation energy needs to be approximated. Approx- 
imations are usually stated by writing either the di- 
agonal of the 2RDM or the xc energy as a functional 
of the natural orbitals and occupation numbers. For 
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most currently employed functionals the xc energy 
takes the form [26, 27, 28, 29, 30, 31] 



1 °° 

9 /K'"'^) 



(34) 



functionals to the time domain leads to a vanishing 
cumulant A. By inserting Eq. (37) into the equation 
of motion for the natural occupation numbers [32], 
we find 



<^j(i-)<^*(r')Vfc(r')¥'fc(r) ,3 ,3 / 

a rd r 



idtUkit) = 'YXijki{t){ij\vint\kl){t) - c.c, 

ijl 



(38) 



r - r' 



i.e. it is given as an exchange integral modified by a 
function depending on the occupation numbers. Gen- 
erally, one can expand the 2RDM in the basis of the 
natural orbitals 

72(ri,r2;r'i,r'2) ^YJ2,^JklVii^^l)Vj{^^2Mi^^'l)^l {'^'2) 

ijkl 

(35) 

which gives rise to expansion coefficients "f2.ijki ■ Re- 
stricting these coefficients to the form 



(36) 



the diagonal, 72(r, r'; r, r'), yields the Hartree energy 
and the xc energy given in Eq. (34). 

Currently, an effort is made to extend the static 
theory in order to describe time-dependent systems. 
The most straightforward extension is again achieved 
by employing an adiabatic approximation [32, 33], i.e. 
the 2RDM at time t is only treated as a functional 
of the IRDM at this point in time. Consequently, 
in Eq. (35), 72 aquires a dependence on time t as do 
the occupation numbers and natural orbitals. For the 
coefficients of the 2RDM (in the basis of the natural 
orbitals at time t) such an adiabatic approximation 
amounts to 

l2,ijki{t) = gijki[{nj{t)}]SikSji ~ hijki[{nj{t)}]5aSjk 
+Kjki[{n,{t)}], (37) 

where gijki (t) and hijki (t) arc the time-dependent co- 
efficients for the Hartree and exchange-type integrals, 
respectively. The time-dependent cumulant Xijki{i) 
contains all contributions to the 2RDM which are not 
of Hartree or exchange type. All three coefficients are 
approximated as functionals of the occupation num- 
bers. Comparing Eqs. (36) and (37) we note, that an 
adiabatic extension of the currently employed static 



which directly illustrates that adiabatic approxima- 
tions based on the form of Eq. (36) cause a zero right- 
hand side and, hence, lead to occupation numbers 
which are constant in time. While one can imagine 
this to be a reasonable approximation in some situa- 
tions, it will generally not be the case. Including an 
explicit cumulant in the approximation of ^2.ijki (t) 
leads to occupation numbers which can aquire a true 
time-dependence, even if one chooses an adiabatic ap- 
proximation for Xijki{t). 

To assess the quality of the adiabatic approxima- 
tion in RDMFT quantitatively, we employ again our 
ID model. For such model systems with a small 
number of electrons one can extract the exact IRDM 
from the solution of the time-dependent Schrodinger 
equation of the interacting system, which allows us 
to explicitly investigate for which situations con- 
stant occupation numbers yield a reasonable descrip- 
tion. For this purpose, we employ two different two- 
electron systems, a ID helium atom and a ID hy- 
drogen molecule. The external potentials are given 

by 



«?xt(2:) 



«cxt(a;) 



(39) 



1 



1 



v/(x + d/2)2 + a2 ^(.x-d/2)2 + a2 

i.e. we are using a soft-Coulomb potential to describe 
the interaction between the nuclei and the electrons 
and, for the hydrogen molecule, also the interaction 
between the two nuclei. 

In the following, we investigate two different 
situations. In the first case we apply a kick, Eq. (9), 
which provides an initial momentum to the system. 
We choose the strength Sq = 0.05 Ha/bohr such that 
the evolution can be described in linear response. As 
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Figure 13: The upper panel displays the laser amplitude of 
an optimized laser pulse that induces a transition from the 
ground state to the first excited singlet state of Helium. In the 
lower panel the two largest occupation numbers of the reduced 
one-body density matrix are shown as function of time. 

a second case, we investigate the transition of the sys- 
tem, here the heUum atom, from its ground state to 
the first excited singlet state. To this end we use op- 
timal control theory [34, 35, 36] to find an optimized 
laser pulse which induces a transition with a popula- 
tion of the excited state of 98.59% at the end of the 
pulse. 

In Fig. 12, we show the dipole moment and the 
change in the first two natural occupation numbers, 
Arij (t) = Uj (t) — Uj (t = 0), for both the helium atom 
and the hydrogen molecule. The strength of the kick 
was chosen as the maximum possible while staying 
within a linear response description. As we can see, 
the occupation numbers show pronounced oscillations 
which, however, remain small in amplitude compared 
to their ground state values while the dipole moment 
shows the characteristic oscillations. Hence, in lin- 
ear response, a description with constant occupation 
numbers will be appropriate. 

As an example where the occupation numbers 



change significantly, we examine the transition of the 
helium atom from its singlet ground state to the first 
excited singlet state which, due to spin, has multi- 
reference character. The optimized laser pulse which 
achieves 98.59% of occupation in the excited state af- 
ter a time of 250 a.u. is shown in Fig. 13. We also 
show the evolution of the first two natural occupation 
numbers which starting close to one and zero, respec- 
tively, approach each other during the propagation. 
In this situation, any description which enforces con- 
stant occupation numbers will clearly not describe 
the situation accurately. 

To describe a situation where occupation numbers 
change in time, functionals in RDMFT have to in- 
corporate time-dependent approximations for the cu- 
mulant of the reduced two-body density matrix [37]. 
Possible approaches along these lines could be based 
on reconstruction approaches for the cumulant of the 
two-body matrix [30], or alternatively on antisym- 
metrized geminal power (AGP) wave functions [38], 
which was proposed recently [32] and will be investi- 
gated in a future study. 

7. Conclusions and outlook 

In this work we aimed to partially unveil the role of 
electron-correlation in the electron dynamics of sys- 
tems driven out of equilibrium. To reach this goal, 
we have used different ID model systems where we 
could assess, by comparing with the exact solution, 
the quality of density and reduced density-matrix 
functionals in various situations where the system 
interacts with external time-dependent fields. We 
looked at both linear and non-linear responses. The 
ID ALDA approximation shows the same behavior 
as its 3D counterpart leading to the well-known un- 
derestimation of the ionization potential and a fail- 
ure to reproduce the excitations to Rydberg states. 
Also, as has been seen in the past [4, 20], we showed 
that in the linear response function double excitations 
do not appear when using adiabatic approximations 
(exemplified here by using both ALDA and EXX). 
For the case studied here, where spatial symmetry 
has been broken and the ground state is a described 
by a doubly occupied KS orbital, double excitations 
become visible beyond linear response for the above 



16 



mentioned adiabatic functionals. In going to the non- 
linear regime, we demonstrate that the description 
of Rabi oscillations within all adiabatic functionals 
leads to a dynamical detuning as the system is driven 
out of resonance by the changes in the potential due 
to the changing density associated with the transition 
during the Rabi oscillation. This manifests itself in 
a very small population of the excited state in con- 
trast to the exact resonant propagation. Hence, the 
description of Rabi oscillations provides a very good 
test case for the development of non-adiabatic func- 
tionals. 

Within RDMFT adiabatic extensions of commonly 
employed ground-state functionals lead to constant 
occupation numbers. This was shown to be a valid 
description within linear response but it turns out to 
be a poor approximation in situations where tran- 
sitions to and among excited states (with possible 
multi-reference character) take place during the evo- 
lution of the system. 

In the future, the ID model systems will be used 
to improve existing approximations, especially go- 
ing beyond the adiabatic dependence on the density 
or the density matrix. Possible routes along these 
lines include e.g. orbital functionals in TDDFT, or 
explicit cumulant approximations in timc-dcpcndcnt 
RDMFT. 
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